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Abstract. We use the Smaller ALignment Index (SALI) method of chaos detection, to study 
the global dynamics of conservative dynamical systems described by differential or difference 
equations. In particular, we consider the well-known 2 and 4-dimensional symplectic standard 
map, as well as an autonomous Hamiltonian system of 2 and 3 degrees of freedom describing 
the motion of a star in models of barred galaxies. The application of SALI helped us to com- 
pute rapidly and accurately the percentage of regular and chaotic motion for particular values 
of the parameters of these systems. We were also able to perform a computationally efficient 
determination of the dependence of these percentages on the variation of several parameters of 
the studied models. 
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^ ■ 1 Introduction 

cn ■ 

The qualitative distinction between chaotic and regular motion in symplectic maps and in systems of 
differential equations is a fundamental problem of non-linear dynamics. This distinction is in general, 
a non trivial task and it becomes more difficult as the number of degrees of freedom increases. For 
this reason, over the years, several methods distinguishing regular from chaotic motion in conservative 
dynamical systems have been proposed and applied, with varying degrees of success. 

One of the most efficient methods of chaos detection is the computation of the so-called Smaller 
ALignment Index (SALI) which was introduced in pQ and has already been applied successfully 
to several dynamical systems [1-18]. SALI has proved to be a fast and reliable method which can 
distinguish between regular and chaotic motion rapidly, reliably and accurately. These characteristics 
^ ■ make this index perfectly suited for the study of global dynamics of dynamical systems. For these 
reasons we use the SALI to study the behavior of two distinct dynamical systems: the well-known 2- 
dimensional (2D) standard map [19] and its generalization to higher dimensions, as well as Hamiltonian 
systems of 2 (2D) and 3 (3D) degrees of freedom, describing the motion of stars in models of barred 
galaxies. In the present paper, we present some preliminary results of our studies. 

The paper is organized as follows: In Section [2] we recall the definition of SALI explaining also 
its behavior for regular and chaotic orbits. In Section [3] we use SALI for computing the fraction of 
chaotic orbits in the case of the standard map, while in Section S] the results of an analogous study for 
Hamiltonian systems of barred galaxies are presented. Finally, in Section[5]we present our conclusions. 



2 The Smaller Alignment Index (SALI) 

Let us consider a I - dimensional phase space of a conservative dynamical system, which could be 
a 2M-dimensional symplectic map or a Hamiltonian flow of N degrees of freedom, with I = 2N. 
We consider also an orbit in this space with initial condition 5(0) = (xx(0), 2:2(0), ...,xi(0)) and two 



deviation vectors Vi(0) = (dxu(0), dxi2(0), dxu(Q)) and V^O) = (dx2i(0), dx22(0), dx2i(0)), from 
the initial point 5(0). 

For the computation of the SALI of a given orbit, one has to follow the time evolution of the orbit 
itself and also of two deviation vectors which initially point in two different directions. The evolution 
of an orbit of a map F is described by the discrete-time equations of the map: 

S(n + l) = F(S(n)), (1) 

where S(n) = (xx(n) , X2(n) , ...,xi(n)), represents the orbit's coordinates at the n-th iteration. The 
evolution of the deviation vectors Vi(n), V2(n), in this case, is given by the equations of the tangent 
map: 

V(n+l) = DF(S(n))-V(n). (2) 
In the case of Hamiltonian flows the evolution of an orbit is given by Hamilton's equations of motion: 

dS(t) 



dt 



F(S(t)), (3) 



where F is a set of n-functions (F\, F2, F n ), while the corresponding evolution of the deviation 
vectors Vi(t),V 2 (t), is given by the variational equations: 

^ = DF(S(t)) ■ V(t) . (4) 

We note that in © and @ DF denotes the Jacobian matrix of equations ([T|) and ^ respectively, 
evaluated at the points of the orbit under study. 

At every time step (or iteration) the two deviation vectors V\{t) and V 2 (t) are normalized with 
norm equal to 1 and the SALI is then computed as: 



SALI it) = min 



Vi(t) + V 2 (t) 



\Vi(t)\\ \\V 2 (t)\\ 



V x {t) V 2 {t) 



\\Vi{t)\\ \\V 2 {t)\\ 



(5) 



where || • || denotes the usual Euclidean norm and t is the continuous or discrete time. 

SALI has a completely different behavior for regular and chaotic orbits and this allows us to clearly 
distinguish between them. In the case of Hamiltonian flows or 2M-dimensional symplectic maps with 
2M > 2, the SALI fluctuates around a non-zero value for regular orbits, while it tends exponentially 
to zero for chaotic orbits [TJH1E], following a rate which depends on the difference of the two largest 
Lyapunov Exponents [6]. Thus, in 2D and 3D Hamiltonian systems the distinction between ordered 
and chaotic motion is easily done. On the other hand, in the case of 2D maps the SALI tends to zero 
both for regular and for chaotic orbits, following however completely different time rates, which again 
allows us to distinguish between the two cases [I] . 



/\ ..;„,.->_ ,. ( m °d 1) , (6) 



3 Global dynamics of 2D and 4D standard map 

As a simple 2D map which exhibits regular and chaotic behavior, we consider the well-known standard 
map [T!5] in the form 

x n +i = x n + y n+ i 

Vn+1 = yn + ^ sin(27TX n ) 

where K is the so-called non-linear parameter of the system. 

Before studying the global dynamics of map ([6]) let us examine in more detail the behavior of SALI 
for regular and chaotic orbits of a 2D map. In the case of a chaotic orbit, any two deviation vectors 
will be aligned to the direction defined by the largest Lyapunov exponent Li, and consequently SALI 
tends to zero following an exponential decay of the form SALI oc e~ 2Lin , with n being the number 
of iterations [6]. In the case of regular orbits any two deviation vectors tend to fall on the tangent 
space of the torus on which the motion lies [H [IJ [18] . For a 2D map this torus is an 1-dimensional 
invariant curve, whose tangent space is also 1-dimensional and consequently any two deviation vectors 
will become aligned. Thus, even in the case of regular orbits in 2D maps the SALI tends to zero. This 
decay follows a power law [T] having the form SALI oc 1/n 2 [18] . 
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Figure 1: The evolution of SALI (solid lines) for a) the chaotic orbit with initial condition xo = 0.2, 
yo = 0.2 and b) the regular orbit with initial condition xo = 0.4, yo = 0.8 of the standard map ([6]) 
for K = 2, with respect to the number of iterations n. Note the different scales of the horizontal 
axis. The Lyapunov exponent of the chaotic orbit is L\ ~ 0.438. Dashed curves in panels a) and b) 
correspond to functions proportional to e~ 2Lin and 1/n 2 respectively. It is evident that the theoretical 
predictions for the evolution of SALI describe very well the numerical data. 

In figure Q] we see the different behavior of SALI for regular and chaotic orbits of the standard 
map d6]). It is exactly this different behavior of the index that allows us to use SALI for a fast and 
clear distinction between regions of chaos and order in the 2-dimensional phase space of the standard 
map. From the results of figure [T] and the theoretical predictions for the evolution of SALI we see that 
after n = 500 iterations the value of SALI of a regular orbit becomes of the order of 10 -6 , while for a 
chaotic orbit SALI has already reached extremely small values. Thus, the percentage of chaotic orbits 
for a given value of K can be computed as follows: We follow the evolution of orbits whose initial 
conditions lie on a 2-dimensional grid of 1000 x 1000 equally spaced points on the 2-dimensional phase 
space of the map (dividing in this way the (x, y) plane in 10 6 squares) and register for each orbit the 
value of SALI after n = 500 iterations. All orbits having values of SALI significantly smaller than 
10 -6 (which correspond to the value SALI reaches after 500 iterations in the case of regular orbits), 
are characterized as chaotic. In practice as a good threshold for this distinction we consider the value 
10 -8 . Thus, all orbits having SALI < 10 -8 after n = 500 iterations are characterized as chaotic, while 
all others are considered as non-chaotic. 

In figure [2^) we present the outcome of this procedure for K = 2. Each initial condition is colored 
according to the color scale seen at the right side of the panel. So, chaotic orbits, having SALI < 10 -8 
are colored black, while light gray color corresponds to regular orbits having high values of SALI. 
Thus, in figure [2^) we can clearly identify even tiny regions of regular motion which are not easily 
seen in phase space portraits of the map (figure Wp))- 

Using the above-described method we were able to compute very fast and accurately the percent- 
ages of regular motion for large values of parameter K. In figured) we plot the percentage of regular 
orbits for 180 < K < 200 where K varies with a step 5K = 0.001. A blow-up of the peak appearing 
close to K = 188 is seen in figure Eb). In order to accelerate the numerical computation we applied 
the following technique: For each orbit we compute its SALI value at n = 500, keeping also track 
of the squares on the (x, y) plane that the orbit visits in its evolution. Then, we attribute the same 
SALI value to all these squares. In this way we gain considerably in computational time, since it is 
not necessary to perform the same computation for the total number of the initial conditions. For 
each value of K a grid of 1000 x 1000 initial conditions were used, allowing us to detect extremely 
tiny regions of regular motion (note that the percentages of regular orbits in figure [3] remain always 
less than 0.0015%!). From the results of figure Owe see a periodicity of period 2tt in the appearance 
of islands of stability as K varies, in accordance to the results presented in [20]. In our study we were 
able to reproduce the results obtained in [20] but with considerably less computational effort. For 
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Figure 2: a) Regions of different values of SALI after n = 500 iterations of map ([75]) for K = 2. b) 
Phase space portrait of map © for the same value of K. 

example, for K = 2 instead of using all the 10 6 initial conditions of the 1000 x 1000 grid, computing 
the evolution of only 12425 initial conditions up to n = 500 iterations was sufficient for characterizing 
the total 10 6 points. Thus, for K = 2 we were able to compute the percentage of regular orbits on 
a 1000 x 1000 grid mesh by computation only 3 x 500 x 12425 « 2 • 10 7 iterations of the map ([6]) 
and its tangent map, instead for the 5 • 10 9 iterations needed for obtaining the same result in |20j . In 
particular for the computation of the data of figure [3] we needed only 27 hours of CPU time on an 
Athlon 64bit, 3.2GHz PC. 

A similar study can be also implemented for the 4D standard map, which is described by the 
following equations: 



where K is again the non-linear parameter and (5 the so-called coupling parameter of the system. In 
figures Hh),b), we present the evolution of the index both for regular and chaotic orbits, for K = 3 
and (3 = 0.1. For the regular orbit of figure 0h.) SALI fluctuates around to a non-zero value, while 
for the chaotic orbit of figure [Da) SALI decays exponentially to zero reaching extremely small values 
after only a few iterations (N ~ 150), following the exponential law: SALI oc e~( Ll ~ L v n : with L\, L2 
being the two largest Lyapunov exponents of the orbit. 

We were also able to measure the percentages of chaotic and regular orbits for the 4D standard 
map following a similar procedure to the one used in the case of the 2D map. We considered 10 6 initial 
conditions equally spaced in the 4-dimensional phase space of the system, producing in this way a fine 
grid of 4-dimensional hypercubes. Noting that in the case of chaotic orbits only a few hundreds of 
iterations are needed for SALI to reach the numerical accuracy of a computer, i.e. SALI ~ 10 -16 (in the 
case of the orbits of figureHb) 150 iterations were sufficient), we started our computation by integrating 
orbits for only 500 iterations. For each orbit we also kept track of the 4-dimensional hypercubes it 
visited in its evolution. If the studied orbit was regular, i.e. its final SALI value was larger than 10~ 8 , 
its final SALI value was attributed to all the hypercubes visited by the orbit. If, on the other hand, 
the orbit was characterized as chaotic (i.e. its SALI value became < 10~ 8 ), the evolution of the orbit 
(but not the evolution of the variational equations) was extended to 5000 iteration, allowing us to 
attribute the computed SALI value to many hypercubes. This procedure decreases significantly the 
CPU time needed for the reliable computation of the percentage of regular motion. In particular for 
K = 3 and (3 = 0.1 the percentage of regular motion was found to be 8,7% after only 1 minute of 
computations with the same computer used in the 2D case. 



x'2 = X2 + j= sin(27rxi) - £ sin[27r(x 3 - sci)] 
x\ = x 4 + J£ sin(27rx 3 ) - ^ sin[27r(xi - x 3 )] 
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Figure 3: a) Percentages of regular orbits of map © as a function of the nonlinear parameter K 6 
[180, 200], b) A zoom of panel a) in the region of K £ (188.44, 188.55). 





Figure 4: The evolution of the SALI for a) the regular orbit with initial condition (x±, X2, X3, X4) = 
(0.55, 0.1, 0.62, 0.2) and b) the chaotic orbit with initial condition (x\, X2, £3, £4) = (0.55, 0.3, 0.62, 0.2) 
of the 4D standard map with K = 3 and (5 = 0.1. We note that the SALI of the chaotic orbit decays 
exponentially to zero following the law oc e~ ( L i- L v n _ 



4 Applications to 2D and 3D models of barred galaxies 

4.1 The model 

A 3D rotating model of a barred galaxy can be described by the Hamiltonian function: 

H =\{pI +pI+pI) + V(x,y,z) - n b (xp y - yp x ). (8) 

The bar rotates around its z-axis, while the x-axis is along the major axis and the y-axis is along 
the intermediate axis. The p x ,Py and p z are the canonically conjugate momenta. Finally, V is the 
potential, U b represents the pattern speed of the bar and H is the total energy of the system. 
The potential V of our model consists of three components: 

1. A disc, represented by a Miyamoto disc |21j : 



+ y 2 + (A + \J z 1 + B 2 ) 2 

where Mjj is the total mass of the disc, A and B are the horizontal and vertical scalelengths, 
and G is the gravitational constant. 



2. A bulge which is modeled by a Plummer sphere whose potential is: 



Vs _ - ■ GM * , (10) 

^Jx 2 + y 2 + z 2 + e 2 s 

where e s is the scalelength of the bulge and Ms is its total mass. 

3. A triaxial Ferrers bar, the density p(x) of which is: 

P ( X ) = { pc{1 - m2)2 ' m<1 , (11) 
PV ' | ,m > 1 V ; 

where p c = IP- is the central density, is the total mass of the bar and 



ijp 1 g 2 

m 2 = -2 + + ^2 > a>6>c>0, (12) 
with a, 6 and c being the semi-axes. The corresponding potential is: 

V B = -nGabc^r H -^-(l - m\u)T + \ (13) 

where 

m2 ^ = J^ u + ¥T^ + ^ (14) 

A 2 (u) = {a 2 + u)(b 2 + u)(c 2 + u), (15) 
n is a positive integer (with n = 2 for our model) and A is the unique positive solution of: 

m 2 {\) = 1, (16) 

outside of the bar (m > 1), and A = inside the bar. 
The corresponding forces are given analytically in |22j . 

This model has been used extensively for orbital studies [23-27]. 
4.2 Numerical results 

We first applied the SALI method to the 2D bar potential resulting from the restriction of our study 
on the z = p z = subspace of the whole phase space of the system. It can be easily seen that, due to 
the symmetries of Hamiltonian ([8]), orbits starting with z = p z = remain for all time on the (x,y) 
plane. In this case, the Hamiltonian function governing the motion is derived by setting z = p z = 
to equations (f8l)- (fT5j) and the corresponding Poincare Surface of Section (PSS) is 2-dimensional and 
can be easily visualized. 

As we have already mentioned, in 2D Hamiltonian systems SALI tends exponentially to zero for 
chaotic orbits, while it fluctuates around a positive number for regular orbits. In figure [5^i) we present 
the PSS (plane (y,p y )) of the system for H = —0.360, which exhibits both regular and chaotic regions. 
By choosing orbits with initial conditions on the line p y = of the PSS and calculating their SALI 
values at t = 3000 we were able to detect very small regions of stability that can not be visualized 
easily on the PSS. We plot the corresponding values of the SALI in figure [5b). The values of the SALI 
tend to zero (~ 10~ 16 ), for orbits whose initial conditions were chosen inside the chaotic regions of the 
PSS, contrary to orbits with initials conditions inside the stability islands whose SALI values retain 
large positive values. 

We also used the SALI method to calculate the percentages of regular and chaotic motion for initial 
conditions chosen on the whole plane (y,p y ). For the value of the Hamiltonian function E = H = 
—0.360 we tested two sets of initial conditions: set A Having 5000 initial conditions and set B with 
10000 initial conditions. The two set were used in order to examine the variation of the percentages 
of the regular and chaotic orbits for different grids of initial conditions. We found that for set A the 




-Oj6 - 



L_i I i I i L_i I i I i_l i I i L_i I i I IU ' ' ' - 

-1.0 -Ofl -OjS -0.4 -02 OJQ 0.2 0.4 OS OS 1.0 08 04 00 

y y 



Figure 5: a) Poincare surface of section for the 2D Ferrers' model in (y,p y ) plane for H = —0.360, b) 
The variation of the SALI value for initial conditions chosen on the line p y = of the corresponding 
PSS of panel a). 

percentage of chaotic orbits is 22, 5%, while for set B 28, 0% of the orbits were characterized as chaotic. 
We note that, as usual, an orbit is considered to be chaotic if its SALI value becomes < 10~ 8 . Thus 
for this particular value of the energy, although a finer grid can help us to detect some small chaotic 
regions, we can nevertheless derive a good approximation of the real percentage even with a relatively 
"small" set of initial conditions. Repeating this procedure for several values of the energy, we were 
able to follow the change of the fraction of chaotic and regular orbits in the phase space as the value 
of H varies. 

We have also studied the complete 3D model, described by the Hamiltonian ([8]), for several values 
of the bar mass Mb and of the semi-minor axis c. Our basic model (main model), has the following 
values of parameters: G = 1, Q b = 0.054, a = 6, b = 1.5, c = 0.6, A = 3, B = 1, e s = 0.4, M B = 0.1, 
M s = 0.08, M D = 0.82 both for its 2D and 3D versions. The units used, are: 1 kpc (length), 1 Myr 
(time), 2x 10 11 solar masses (mass). For each model studied we considered two sets of initial conditions. 
The first set contains orbits with initial conditions in the (x,p y ,z) space with (y,p x ,Pz) = (0,0,0) 
and the second one, orbits with initial conditions in the (x,p y ,p z ) space with (y,z,p x ) = (0,0,0). 
Analyzing our numerical results we found that the increase of the bar mass (2GMB - version, with 
Mb = 0.2) introduces more chaotic behavior for both sets of initial conditions. In figure [6] we present 
the change of percentages of the chaotic orbits as the parameters of Hamiltonian (|8|) vary for the 
first set of initial conditions (similar results where also found for the second set). The findings are 
in accordance to the results obtained in [28J for the 2D case. On the other hand, we discovered that 
when the bar is thicker, i.e. the length of the z-axis larger (2C - version, with c = 1.2), the system 
becomes less chaotic [TT] . 

Finally, we calculated the percentages of chaotic and regular orbits for different values of the 
pattern speed f2&. From the orientation of periodic orbits, Contopoulos [29] showed that bars have to 
end before corotation, i.e. that > a, where tl the Lagrangian, or corotation, radius. Comparing 
the shape of the observed dust lanes along the leading edges of bars to that of the shock loci in 
hydrodynamic simulations of gas flow in barred galaxy potentials, Athanassoula [30] [31] was able to 
set both a lower and an upper limit to corotation radius, namely tl = (1-2 ± 0.2)a. This restricts the 
range of possible values of the pattern speed between a high value that corresponds to the Lagrangian 
radius tl = 1.4a and a low value that corresponds to ri, = 1.0a. Using the extremes of this range, 
we investigated how the pattern speed of the bar affects the dynamics the system and found that the 
percentage of regular orbits is greater in slow bars |16j . 




Figure 6: Comparison of the percentages of chaotic orbits for our main galactic model, described by 
the Hamiltonian function ([8]) with parameters: G = 1, f2j, = 0.054, a = 6, b = 1.5, c = 0.6, ^4 = 3, 
B = 1, e s = 0.4, M B = 0.1, M s = 0.08, M D = 0.82, a model with twice the length of the short axis 
c = 1.2 (#C - version) and a model with twice the bar mass Mb = 0.2 (2GMB - version). The system 
becomes more chaotic as the mass of the bar component increases, while a thicker bar results to the 
decrease of chaoticity. The initial conditions in this example are given in the (x,p y ,z) space with 
{y,Px,Pz) = (0,0,0) 



5 Conclusions 

In this paper, we used the SALI method of chaos detection to study the dynamical behavior of 2 and 
4 dimensional symplectic maps and of Hamiltonian models of barred galaxies with 2 and 3 degrees of 
freedom. Using the SALI we were able to rapidly identify, in the phase spaces of the studied systems, 
even tiny regions of regular motion and also compute the percentages of regular and chaotic orbits as 
the values of the parameters of the systems vary. In the case of galactic models, we found the influence 
of some important physical parameters like the mass, the length of the short z-axis and the pattern 
speed of the bar, on the chaotic behavior of the system. In particular, we found that the growth of 
the mass of the bar favors the existence of more chaotic orbits, while we observed that by increasing 
the length of the short axis of the bar the percentage of the chaotic orbits decreases. 
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